%dynaremodel_KurmannSims_Part3.mod
%
%dynare file called to solve and simulate the NK DSGE model

var G Gh Ge lam1 lam2 vap C int Pi R z theta e h mrs mui I K Lam f1 f2 w Ld N wr x1 x2 pm Y pr A Kser vp L vw g gX du dh duob dtfp dtfpu dtfpc dY dC dK dhN hN hu dtech lshare;

varexo eA ei eg et em ep eu;

parameters Gs alpha beta delta phiw phip gammaw gammap epsp epsw gXs b Pis k0 k1 k2 k3 k4 gamma1 gamma2 phi psi thetas T Ns iks gs rhoi phipi phiy si rhoA rhog rhot rhom rhop sA sg st sm sp hs es lam1s Cs Ys vps vws Is Ks ints Rs ws wrs mrss prs pms Lams betah su rhou F;

load param_new_fixed_new;

set_param_value('Gs',Gs);
set_param_value('alpha',alpha);
set_param_value('beta',beta);
set_param_value('delta',delta);
set_param_value('phiw',phiw);
set_param_value('phip',phip);
set_param_value('gammaw',gammaw);
set_param_value('gammap',gammap);
set_param_value('epsp',epsp);
set_param_value('epsw',epsw);
set_param_value('gXs',gXs);
set_param_value('b',b);
set_param_value('Pis',Pis);
set_param_value('k0',k0);
set_param_value('k1',k1);
set_param_value('k2',k2);
set_param_value('k3',k3);
set_param_value('k4',k4);
set_param_value('gamma1',gamma1);
set_param_value('gamma2',gamma2);
set_param_value('phi',phi);
set_param_value('psi',psi);
set_param_value('thetas',thetas);
set_param_value('T',T);
set_param_value('Ns',Ns);
set_param_value('iks',iks);
set_param_value('gs',gs);
set_param_value('rhoi',rhoi);
set_param_value('phipi',phipi);
set_param_value('si',si);
set_param_value('rhoA',rhoA);
set_param_value('rhog',rhog);
set_param_value('rhot',rhot);
set_param_value('rhom',rhom);
set_param_value('rhop',rhop);
set_param_value('sA',sA);
set_param_value('sg',sg);
set_param_value('sp',sp);
set_param_value('sm',sm);
set_param_value('st',st);
set_param_value('phiy',phiy);
set_param_value('hs',hs);
set_param_value('es',es);
set_param_value('lam1s',lam1s);
set_param_value('Ys',Ys);
set_param_value('Cs',Cs);
set_param_value('Is',Is);
set_param_value('Ks',Ks);
set_param_value('Rs',Rs);
set_param_value('ints',ints);
set_param_value('vps',vps);
set_param_value('vws',vws);
set_param_value('ws',ws);
set_param_value('wrs',wrs);
set_param_value('mrss',mrss);
set_param_value('prs',prs);
set_param_value('pms',pms);
set_param_value('Lams',Lams);
set_param_value('betah',betah);
set_param_value('su',su);
set_param_value('rhou',rhou);
set_param_value('F',F);

model;
% (1) G function
G = k0 + exp(hu)*(k1/k2)*exp(h)^(k2) + (k3/k4)*exp(e)^(k4);                  

% (2) Gh (partial wrt h)
Gh = exp(hu)*k1*exp(h)^(k2-1);

% (3) Ge  (partial wrt e)
Ge = k3*exp(e)^(k4-1);

% (4) lam1 (MUC)
exp(lam1) = exp(vap)/(exp(C) - b*gX^(-1)*exp(C(-1))) - beta*b*exp(vap(+1))/(gX(+1)*exp(C(+1)) - b*exp(C));

% (5) Euler equation bonds
exp(lam1) = beta*(1+int)*gX(+1)^(-1)*exp(lam1(+1))*exp(Pi(+1))^(-1);

% (6) Utilization FOC
exp(R) = gamma1 + gamma2*(exp(z) - 1);

% (7) FOC e 
exp(theta)*exp(vap)*Ge = exp(lam1)*exp(mrs)*exp(h);

% (8) FOC h
exp(theta)*exp(vap)*Gh = exp(lam1)*exp(mrs)*exp(e);

% (9) Investment FOC
exp(lam1) = exp(lam2)*exp(mui) - phi*(exp(I)/exp(K(-1))*gX - iks)*exp(lam1);

% (10) N foc
exp(theta)*exp(vap)*G = exp(lam1)*(exp(mrs)*exp(e)*exp(h) - (psi/2)*(exp(N)/exp(N(-1)) - 1)^2*exp(mrs) - psi*(exp(N)/exp(N(-1)) - 1)*exp(mrs)*exp(N)/exp(N(-1))) + beta*exp(lam1(+1))*exp(mrs(+1))*psi*(exp(N(+1))/exp(N) - 1)*(exp(N(+1))/exp(N))^2;

% (11) K' FOC 
exp(lam2) = beta*gX(+1)^(-1)*(exp(lam1(+1))*(exp(R(+1))*exp(z(+1)) - gamma1*(exp(z(+1)) - 1) - (gamma2/2)*(exp(z(+1)) - 1)^2 - (phi/2)*(exp(I(+1))/exp(K)*gX(+1) - iks) + (exp(I(+1))/exp(K)*gX(+1))*phi*(exp(I(+1))/exp(K)*gX(+1) - iks)) + exp(lam2(+1))*(1-delta));

% (12) SDF
Lam = beta*gX^(-1)*exp(lam1)/exp(lam1(-1));

% (13) f1 (wage)
exp(f1) = exp(mrs)*exp(w)^(epsw)*exp(Ld) + phiw*gXs^(-epsw - 1)*gX(+1)^(1+epsw)*Lam(+1)*exp(Pi)^(-epsw*gammaw)*exp(Pi(+1))^(epsw)*exp(f1(+1));

% (14) f2 (wage)
%exp(f2) = exp(w)^(epsw)*exp(Ld) + phiw*gX(+1)^(epsw)*Lam(+1)*exp(Pi)^((1-epsw)*gammaw)*exp(Pi(+1))^(epsw-1)*exp(f2(+1));
exp(f2) = exp(w)^(epsw)*exp(Ld) + phiw*gXs^(-epsw)*gX(+1)^(epsw)*Lam(+1)*exp(Pi)^((1-epsw)*gammaw)*exp(Pi(+1))^(epsw-1)*exp(f2(+1));


% (15) wage reset
wr = (epsw/(epsw-1))*exp(f1)/exp(f2);

% (16) x1 (price)
x1 = exp(pm)*exp(Y) + phip*gX(+1)*Lam(+1)*exp(Pi)^(-epsp*gammap)*exp(Pi(+1))^(epsp)*x1(+1);

% (17) x2 (price)
x2 = exp(Y) + phip*gX(+1)*Lam(+1)*exp(Pi)^((1-epsp)*gammap)*exp(Pi(+1))^(epsp-1)*x2(+1);

% (18) Reset price
pr = (epsp/(epsp-1))*x1/x2;

% (19) Production
exp(Y) = exp(A)*exp(Kser)^(alpha)*exp(Ld)^(1-alpha) - F;

% (20) FOC capital
exp(R) = alpha*exp(pm)*exp(A)*exp(Kser)^(alpha-1)*exp(Ld)^(1-alpha);

% (21) FOC labor
exp(w) = (1-alpha)*exp(pm)*exp(A)*exp(Kser)^(alpha)*exp(Ld)^(-alpha);

% (22) Output and wholesale output
%exp(Ym) = exp(Y)*vp;

% (23) Price dispersion
vp = (1-phip)*pr^(-epsp) + phip*exp(Pi(-1))^(-gammap*epsp)*exp(Pi)^(epsp)*vp(-1);

% (24) Price evolution
1 = (1-phip)*pr^(1-epsp) + phip*exp(Pi(-1))^(gammap*(1-epsp))*exp(Pi)^(epsp-1);

% (25) Labor market clearing
exp(L) = exp(Ld)*vw;

% (26) Wage dispersion
vw = (1-phiw)*(wr/exp(w))^(-epsw) + phiw*gX(+1)^(epsw)*gXs^(-epsw)*exp(Pi(-1))^(-gammaw*epsw)*exp(Pi)^(epsw)*(exp(w)/exp(w(-1)))^(epsw)*vw(-1);

% (27) Wage evolution
exp(w)^(1-epsw) = (1-phiw)*wr^(1-epsw) + phiw*gX(+1)^(epsw-1)*gXs^(1-epsw)*exp(Pi(-1))^(gammaw*(1-epsw))*exp(Pi)^(epsw-1)*exp(w(-1))^(1-epsw);

% (28) Labor supply
exp(L) = exp(e)*exp(h)*exp(N);

% (29) Capital services and capital
exp(Kser) = exp(z)*exp(K(-1))*gX^(-1);

% (30) Resource constraint
exp(Y) = exp(C) + exp(I) + (psi/2)*(exp(N)/exp(N(-1)) - 1)^2*exp(mrs)*exp(N) + (phi/2)*(exp(I)/exp(K(-1))*gX - iks)^2*exp(K(-1))*gX^(-1) + (gamma1*(exp(z)-1) + (gamma2/2)*(exp(z)-1)^2)*exp(K(-1))*exp(gX)^(-1);

% (31) Capital accumulation
exp(K) = exp(mui)*exp(I) + (1-delta)*exp(K(-1))*gX^(-1);

% (32) Taylor rule
int = (1-rhoi)*ints + rhoi*int(-1) + (1-rhoi)*(phipi*(Pi - log(Pis)) + phiy*(Y - Y(-1) + log(gX) - log(gXs))) + si*ei;

% (33) Process for g  
g = (1-rhog)*log(gs) + rhog*g(-1) + sg*eg(-1);
%g = (1-rhog)*log(gs) + rhog*g(-1) + sg*eg;         // USE THIS IF NEWS SHOCK HAS CONTEMPORANEOUS EFFECT

% (34) gX
gX = exp(g)^(1/(1-alpha));

% (35) Process for A
A = rhoA*A(-1) + sA*eA;

% (36) Process for theta
theta = (1-rhot)*log(thetas) + rhot*theta(-1) + st*et;

% (37) Process for mui
mui = rhom*mui(-1) + sm*em;

% (38) Process for vap
vap = rhop*vap(-1) + sp*ep;

% (39) du (change in actual utilization)
du = alpha*(z - z(-1)) + (1-alpha)*(e - e(-1));

% (40) dh (change in hours)
dh = h - h(-1);

% (41) duob (change in observed utilization)
duob = betah*dh;

% (42) Labor share
lshare = (exp(w)*exp(Ld))/exp(Y);

% (43) dtfp
dtfp = dY - (1-lshare)*dK(-1) - lshare*dhN;

% (44) dtfpu
dtfpu = dtfp - du;

% (45) dtfpc
dtfpc = dtfp - duob;

% (46) dY
dY = Y - Y(-1) + log(gX);

% (47) dC
dC = C - C(-1) + log(gX);

% (48) dK
dK = K - K(-1) + log(gX);

% (49) dhN
dhN = h + N - (h(-1) + N(-1));

% (50) hn
hN = h + N;

% (51) hu (h labor supply shifter)
hu = rhou*hu(-1) + su*eu;

% (52) dtech
dtech = g + A - A(-1);

end;

initval;
G = Gs;
Gh = k1*hs^(k2-1);     
Ge = k3*es^(k4-1);     
lam1 = log(lam1s);
lam2 = log(lam1s);
C = log(Cs);
Y = log(Ys);
e = log(es);
h = log(hs);
N = log(Ns);
I = log(Is);
K = log(Ks);
Kser = log(Kss);
int = ints;
R = log(Rs);
vw = vws;
vp = vps;
L = log(es*hs*Ns);
Ld = log(es*hs*Ns/vws);
w = log(ws);
wr = wrs;
mrs = log(mrss);
pr = prs;
pm = log(pms);
x1 = pms*Ys/(1-phip*gXs*Lams*Pis^(epsp*(1-gammap)));
x2 = Ys/(1-phip*gXs*Lams*Pis^((1-epsp)*(gammap-1)));
f1 = log(mrss*ws^(epsw)/(1-phiw*Lams*Pis^(epsw*(1-gammaw))));
f2 = log(ws^(epsw)/(1-phiw*Lams*Pis^((1-epsw)*(gammaw-1))));
A = 0;
vap = 0;
mui = 0;
theta = log(thetas);
g = log(gs);
z = 0;
Lam = Lams;
gX = gXs;
Pi = log(1);
du = 0;
duob = 0;
dh = 0;
dtfp = log(gs);
dtfpc = log(gs);
hN = log(hs) + log(Ns);
lshare = ws*Lds/Ys;
end;


%steady(solve_algo=0);
steady;

shocks;
var eA = 1;
var ei = 1;
var et = 1;
var ep = 1;
var em = 1;
var eg = 1;
var eu = 1;
end;

stoch_simul(order=1,irf=200,nograph,ar=0);